home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / AMIGA / AMICUS / AMICUS18.ADF / Progs / StarProbe / SPINIT.C < prev    next >
Text File  |  1989-01-27  |  6KB  |  241 lines
  1. /****************************************************************************
  2.  ***                                                                      ***
  3.  ***      S T A R   P R O B E       I N I T I A L I Z E                   ***
  4.  ***                                                                      ***
  5.  ****************************************************************************/
  6. /*****
  7.  ***  This module opens libraries, opens files, gets options, and
  8.  ***  gets initial parameters.
  9.  *****/
  10.  
  11. #include "stdio.h"
  12. #include "math.h"
  13. #include "limits.h"
  14. #include <dos.h>
  15. #include "spmac.h"
  16. #include "spref.h"
  17.  
  18. extern FILE *pf;
  19.  
  20. /*****
  21.  *** I n i t i a l i z e
  22.  *****/
  23.  
  24. void initialize()
  25. {
  26. iam(1);
  27. int i,j;
  28. FILE *inf;
  29. double d;
  30.  
  31.  pf = fopen("sprpt.rpt","w");
  32.  inf= fopen("spparm.in","r");
  33.  
  34.  printf("**************************************************\n");
  35.  printf("***        S  T  A  R  P  R  O  B  E           ***\n");
  36.  printf("***                                            ***\n");
  37.  printf("***        version 1.0    19-nov-86            ***\n");
  38.  printf("**************************************************\n");
  39.  fprintf(pf,"**************************************************\n");
  40.  fprintf(pf,"***        S  T  A  R  P  R  O  B  E           ***\n");
  41.  fprintf(pf,"***                                            ***\n");
  42.  fprintf(pf,"***        version 1.0    19-nov-86            ***\n");
  43.  fprintf(pf,"**************************************************\n");
  44.  
  45.  
  46. /* set defaults */
  47.  
  48.  for (i=0;i<MAX_PROCS;i++) proc_trace[i]=0;
  49.  massratio = 1.0;
  50.  debug = 0;
  51.  max_tries = 2;
  52.  X = 0.70;  Y = 0.27;  Z = 0.03;
  53.  fitptmass = 0.5;
  54.  massgrad[ 0] = 0.001;
  55.  massgrad[ 1] = 0.001;
  56.  massgrad[ 2] = 0.001;
  57.  massgrad[ 3] = 0.001;
  58.  massgrad[ 4] = 0.001;
  59.  massgrad[ 5] = 0.001;
  60.  massgrad[ 6] = 0.001;
  61.  massgrad[ 7] = 0.001;
  62.  massgrad[ 8] = 0.001;
  63.  massgrad[ 9] = 0.001;
  64.  massgrad[10] = 0.01;
  65.  massgrad[11] = 0.01;
  66.  massgrad[12] = 0.01;
  67.  massgrad[13] = 0.01;
  68.  massgrad[14] = 0.01;
  69.  massgrad[15] = 0.01;
  70.  massgrad[16] = 0.01;
  71.  massgrad[17] = 0.01;
  72.  massgrad[18] = 0.01;
  73.  massgrad[19] = 0.05;
  74.  massgrad[20] = 0.05;
  75.  massgrad[21] = 0.10;
  76.  massgrad[22] = 0.10;
  77.  massgrad[23] = 0.10;
  78.  massgrad[24] = 0.10;
  79.  massgrad[25] = 0.10;
  80.  massgrad[26] = 0.10;
  81.  massgrad[27] = 0.05;
  82.  massgrad[28] = 0.05;
  83.  massgrad[29] = 0.01;
  84.  massgrad[30] = 0.01;
  85.  massgrad[31] = 0.01;
  86.  massgrad[32] = 0.01;
  87.  massgrad[33] = 0.01;
  88.  massgrad[34] = 0.01;
  89.  massgrad[35] = 0.01;
  90.  massgrad[36] = 0.01;
  91.  massgrad[37] = 0.01;
  92.  massgrad[38] = 0.01;
  93. /*
  94.  massgrad[38] = 0.001;
  95.  massgrad[39] = 0.001;
  96.  massgrad[40] = 0.001;
  97.  massgrad[41] = 0.001;
  98.  massgrad[42] = 0.001;
  99.  massgrad[43] = 0.001;
  100.  massgrad[44] = 0.001;
  101.  massgrad[45] = 0.001;
  102.  massgrad[46] = 0.001;
  103.  massgrad[47] = 0.001;
  104. */
  105.  gradsteps = 39;
  106.  
  107. /* get parmfile data */
  108.  
  109.  if (fscanf(inf,"%lf ",&massratio)==EOF) goto EndInit;
  110.  
  111.  if (fscanf(inf,"%d ",&debug)==EOF) goto EndInit;
  112.  
  113.  if (debug==1){
  114.    i = 99;
  115.    while (i > 0){
  116.      if (fscanf(inf,"%d ",&i)==EOF) goto EndInit;
  117.      if ((i > 0) & (i < MAX_PROCS)) proc_trace[i] = 1;
  118.      }
  119.    }
  120.  
  121.  if (debug>0) printf("Debug = %d\n",debug);
  122.  
  123.  block(in,"Initialize",0.0);
  124.  
  125.  if (fscanf(inf,"%d ",&max_tries)==EOF) goto EndInit;
  126.  
  127.  if (fscanf(inf,"%lf ",&X)==EOF) goto EndInit;
  128.  if (fscanf(inf,"%lf ",&Y)==EOF) goto EndInit;
  129.  if (fscanf(inf,"%lf ",&Z)==EOF) goto EndInit;
  130.  
  131.  if (fscanf(inf,"%lf ",&fitptmass)==EOF) goto EndInit;
  132.  if (fscanf(inf,"%d ",&gradsteps)==EOF) goto EndInit;
  133.  d=0.0;
  134.  for (i=0;i<gradsteps;i++)
  135.    if (fscanf(inf,"%lf ",&massgrad[i])==EOF){
  136.      printf("*** Insufficient mass elements = %d",i);
  137.      exit(1);
  138.      }
  139.    else d+=massgrad[i];
  140.  
  141.  
  142. EndInit:
  143.  mass=SM*massratio;
  144.  u=meanmolweight(X,Y,Z);
  145.  
  146.  printf("Massratio = %g\n",massratio);
  147.  printf("Trials    = %d\n",max_tries);
  148.  printf("Gradsteps = %d\n",gradsteps);
  149.  printf("FitPtMass = %g\n",fitptmass);
  150.  printf("Total mass fraction = %g\n",d);
  151.  printf("X=%g, Y=%g, Z=%g\n",X,Y,Z);
  152.  if (debug > 0) fprintf(pf,"\n Debug trace on. Level=%d.\n",debug);
  153.  
  154.  for (i=0;i<4;i++)
  155.    for (j=0;j<6;j++)
  156.      results[i][j] = 0.0;
  157.  
  158.  adj_p = 1.001;
  159.  adj_t = 1.001;
  160.  adj_l = 1.001;
  161.  adj_r = 1.001;
  162.  cur_diff = 0.0;
  163.  old_diff = 0.0;
  164.  
  165.  fclose(inf);
  166.  
  167.  block(out,"Initialize",mass);
  168. }
  169.  
  170. /*****
  171.  *** S e t   B o u n d a r y   C o n d i t i o n s
  172.  ***
  173.  *** Ref: Bowers Eq. 9.9 - 9.23 
  174.  *****/
  175.  
  176. void set_boundary_conditions()
  177. {
  178.  iam(2);
  179.  block(in,"set_boundary_conditions",massratio);
  180.  
  181.  effpress = 6.0e9;
  182.  if (massratio == 1.0){ /* special test case */
  183.    radratio = 0.94;
  184.    lumratio = 0.72;
  185.    coredensity = 90.0;
  186.    coretemp = 13.7e6;
  187.    efftemp = 5900.0;
  188.    effpress = 6.0e9;
  189.    }
  190.  else if (massratio<=2.0){/*p-p cycle*/
  191.    radratio=0.312*pow(u,-0.538)*pow(massratio,0.0769);
  192.    lumratio=49.1*pow(u,7.77)*pow(massratio,5.46);
  193.    coretemp=3.05e7 * pow(u,1.54)*pow(massratio,0.923);
  194.    coredensity=86.0*pow(u,-1.615)*pow(massratio,0.769);
  195.    efftemp=2.73e4 * pow(u,2.21)*pow(massratio,1.327);
  196.    }
  197.  else if (massratio<=3.0){/*CNO cycle*/
  198.    radratio=0.451*pow(u,0.395)*pow(massratio,0.697);
  199.    lumratio=43.5*pow(u,7.30)*pow(massratio,5.18);
  200.    coretemp=1.98e7 * pow(u,0.606)*pow(massratio,0.364);
  201.    coredensity=65.8*pow(u,-0.455)*pow(massratio,-0.909);
  202.    efftemp=2.34e4 * pow(u,1.63)*pow(massratio,0.871);
  203.    }
  204.  else{/*electron scattering*/
  205.    radratio=0.454*pow(u,0.588)*pow(massratio,0.765);
  206.    lumratio=112.0*pow(u,4.0)*pow(massratio,3.0);
  207.    coretemp=2.12e7 * pow(u,0.412)*pow(massratio,0.235);
  208.    coredensity=60.3*pow(u,-1.765)*pow(massratio,-1.294);
  209.    efftemp=2.77e4 * pow(u,0.706)*pow(massratio,0.368);
  210.    }
  211.    /* pressure from Clayton eq 2.10 */
  212.  corepress=N0*K*coredensity*coretemp/u;
  213.  B = 0.8;      /*ratio of gas pressure to total pressure */
  214.  poly = 1.25;  /* polytropic index */
  215.  
  216.  block(out,"set_boundary_conditions",u);
  217. }
  218.  
  219.  
  220. /*****
  221.  *** S h u t d o w n
  222.  *****/
  223.  
  224. void shutdown()
  225. {
  226. iam(99);
  227. int i;
  228.  block(in,"shutdown",0.0);
  229.  
  230.  report_results();
  231.  
  232.  for (i=0;i<MAX_PROCS;i++)
  233.    if (proc_count[i] > 0)
  234.     fprintf(pf,"proc %d = %d\n",i,proc_count[i]);
  235.  
  236.  block(out,"shutdown",0.0);
  237.  
  238.  fclose(pf);
  239. }
  240.  
  241.